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An open question is whether the liquid and glassy phases of water are ther- 
mo dynamically distinct or continuous Here we address this question using 
molecular dynamics simulations in comparison with neutron scattering experi- 
ments to study the effect of temperature and pressure on the local structure 
of liquid water. From both simulations and experiments, we find that the 
liquid structure at high pressure is nearly independent of temperature, and 
remarkably similar to the known structure of the high-density amorphous ice 
(HDA). Further at low pressure, the liquid structure appears to approach the 
experimentally- measured structure of low-density amorphous ice (LDA)as tem- 
perature decreases. These results are consistent with continuity between the 
liquid and glassy phases of H 2 0. 

The structure of liquid water has been well-studied at ambient pressure by a variety of 
experimental and simulation techniques. It has been recognized that each water molecule 
is typically hydrogen bonded to four neighboring molecules in a tetrahedral arrangement, 
leading to an open bond network that can account for a variety of the known anomalies of 
water More recently, the effect of pressure on both the structure and the hydrogen bond 
network of liquid water has been studied experimentally @-§j] and by simulations using a 
variety of potentials, including the ST2 potential the MCY potential [[|, the TIP4P 

potential |]7|j8|J10|--|T2| , and the SPC/E potential [pT3|,p~4|| . Furthermore, understanding the 
effects of pressure may be useful in elucidating the puzzling behavior of liquid water. 



In particular, three competing "scenarios" have been hypothesized to explain the anoma- 
lous properties of water: (i) the existence of a spinodal bounding the stability of the liquid 
in the superheated, stretched, and supercooled states |L^,[nj; (ii) the existence of a liquid- 
liquid transition line between two liquid phases differing in density [P^ fr^fi7|fr8[| ; (iii) a 



singularity-free scenario in which the anomalies are related to the presence of low-density 
and low-entropy structural heterogeneities fl9fl . 

Here, we present molecular dynamics (MD) simulations (Table p of a comparatively 
large system of 8000 molecules interacting via the extended simple point charge (SPC/E) 
pair potential [p4j] . We find remarkable agreement with neutron scattering studies of the 



effect of pressure on the structure of liquid D 2 [§]], indicating that the SPC/E potential 
reproduces many structural changes in the liquid over a wide range of temperature and 
pressure. By comparing the simulated pair correlation functions and structure factor with 
our experimental data, we find that the structure of the supercooled liquid at low pressure 
resembles the structure of low-density amorphous (LDA) ice. At high pressure, we find that 
the structure of the liquid appears independent of temperature and is nearly indistinguish- 
able from that of high-density amorphous (HDA) ice. The combined results at high and low 
pressure appear consistent with continuity between the liquid and glassy states of water and 
also with a possible liquid-liquid transition in the supercooled region of the liquid terminated 
by a second critical point. 

We analyze the structures found in our simulations by calculating the atomic radial 
distribution functions (RDF) goo( r ), 9oH( r ), and guiiix)- To compare the RDFs with 
neutron scattering measurements, we form the weighted sum 

h(r) = A-npr [wig o(r) + w 2 goH{r) + w 3 g H H(r) - 1] , (1) 

where the weighting factors Wi are selected to coincide with experimental measurements of 
D2O. Experimentally, h{r) is obtained by Fourier transformation of the molecular structure 
factor. The dominant contributions to h(r) are the H-H and O-H (or D-D and O-D for 
D2O) spatial correlations, so h(r) includes relatively little contribution from oxygen-oxygen 



correlations. Fig. |l|a shows h(r) at two of the five temperatures simulated and also compares 
with experimental data. The peaks centered at 1.8 A, 2.3 A, and 2.8 A correspond to the 
O-H, H-H, and 0-0 intermolecular distances in the hydrogen-bonded configuration, respec- 
tively. While the magnitudes of these peaks change slightly, their ubiquity demonstrates the 
stability of the first neighbor ordering - namely that each molecule is typically surrounded 
by four molecules in a tetrahedral configuration. 

We find that the peak at 3.3 A of h(r) becomes more pronounced as P increases. Ex- 
amination of the individual RDFs shows that the increase at 3.3 A can be attributed to 
changes in garniy). In addition, goo( r ) shows a pronounced increase at 3.3 A under pressure 
(Fig. |l]b), but cannot account for the changes in h(r), as the goo( r ) weighting factor in 
Eq. (JJ) is small [0-0]. The growth at 3.3 A in goo( r ) indicates that the liquid locally has 
the structure of an interpenetrating tetrahedral network similar to ice VI and VII, the ice 
polymorphs close to the high pressure liquid, and can also be associated with the formation 
of clusters with structure similar to HDA [^5| . 

To directly compare with experimental measurements, we calculate the molecular struc- 
ture factor S M (q) = S % ^ tra (q) + S l ^ ter (q), where S l ^ tra (q) and S t ^j ter (q) are the intramolecular 
and intermolecular contributions to Sm{q)i respectively. We calculate S l ^} tra (q) explicitly, 



as described in ref. ||26|| , and Fourier transform h(r) to obtain S^} ter (q). We find striking 
agreement between simulated and the experimental values of SM{q) (Fig. |I]c). In particular, 
the value q = q (P,T) of the first peak of Sm(q) shifts to larger q values as P increases 
(Fig. 0). We note the feature that as temperature is lowered: (i) for P = 600 MPa, q 
approaches the experimental go value of HDA ice (2.20 A -1 ), (ii) For P = 0.1 MPa and 
—200 MPa, qo approaches the experimental qo value of LDA ice (1.69 A -1 ). Indeed we find 
the high-pressure liquid structure resembles that of HDA ice (Fig. |3|a), and the low-pressure 
liquid structure that of LDA ice (Fig. |3|b) . 

We also study the structure of the glassy phase and find that the structure of the sim- 
ulated glasses strongly resembles the experimentally-measured structure of HDA and LDA 
solid water (Fig. H, c and d), suggesting that the simulated glassy state points are the 



analogues of HDA and LDA solid water for SPC/E. Furthermore, by compression of the 
simulated LDA-like system, we are able to reversibly transform the structure to the HDA- 
like system, reminiscent of the experimentally observed reversible first-order transition of 
LDA to HDA under pressure p9| . 

These results support the postulated continuity between atmospheric pressure water and 
LDA ice, based on measurements of the free energy |K| and dielectric relaxation time at 
1 atm. Our results are also consistent with continuity between high pressure liquid water and 
HDA ice [HHH. Combined with the experimentally-detected first-order transition between 
HDA and LDA, our results are consistent with the presence the postulated second critical 
point in the supercooled region of the phase diagram. 



Methods 

We equilibrate systems of 8000 molecules interacting via the SPC/E potential to a constant tem- 
perature and pressure by monitoring the evolution of the density and internal energy. The SPC/E 
model treats water as a rigid molecule consisting of three point charges located at the atomic cen- 
ters of the oxygen and hydrogen which have an OH distance of 1.0 A and HOH angle of 109.47°, 
the tetrahedral angle. Each hydrogen carries a charge qu = 0.4238e and the oxygen carries a 
charge qo = —2qH, where e is the magnitude of the electron charge. In addition, the oxygen atoms 
of separate molecules interact via a Lennard-Jones potential with parameters a = 3.166 A and 



e = 0.6502 kJ/mol. We adjust the temperature and pressure via the methods of Berendsen [2C] 
and we account for the long-range Coulomb interactions using the reaction field technique with a 
cutoff of 0.79 nm pi] . The equations of motion evolve by the SHAKE algorithm |2^| with a time 
step of 1 fs. We typically simulate each state point using eight processors in parallel. We obtain 
simulation speeds of approximately 15/is per particle per update. The total simulation time is 
about 1000 CPU-days. 
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TABLES 

TABLE I. Summary of results from the 24 liquid and 2 glassy state points simulated. The 
liquid state points simulated at positive pressure correspond roughly to temperatures and pres- 
sures studied in the experiments of ref. Q. To facilitate comparison between simulations and 
experiments, we define AT = T — Tmd, the temperature relative to that of the 1 atm tempera- 
ture of maximum density Tmd- At atmospheric pressure, the SPC/E potential displays a density 
maximum at T^fS^ E ~ 245 K fl3],14,2"3||. Similarly, experimental temperatures reported relative 



to the 1 atm Tmd of D2O, 284 K. All state points are liquid, except AT = —145 °C, which is 
glassy. We obtain state points at AT = —145 °C and P = 0.1 and 600 MPa by quenching a 
configuration of the supercooled liquid state points at AT = —30 °C. The configurations ana- 
lyzed at AT = —145 °C are not "equilibrated", as the required simulation time far exceeds the 
computational resources available; rather, they correspond to a glassy state quenched from the 
supercooled liquid. We consider these state points to be glassy since previous simulations indicate 



that SPC/E approaches a glass for AT < —50 °C [23,2q|. At negative pressures, smaller systems 
may not reproduce cavitation events that we observe in larger systems. Thus it is important that 
we consider the large 8000 molecule systems for P = —200 MPa. We also simulated one system of 
64,000 molecules, but did not observe any significant differences in structure or cavitation events. 



P (MPa) 


AT (°C) 


p (g/cm 3 ) 


U (kJ/mol) 


Equilibration Time (ps) 


600 


35 


1.1930 ± 0.0003 


-47.65 


1100 




10 


1.2041 ± 0.0004 


-49.02 


1300 




-10 


1.2139 ±0.0005 


-50.22 


2000 




-30 


1.2236 ± 0.0003 


-51.43 


3000 




-145 


1.250 ±0.008 


-55.3 


4000 


465 


35 


1.1616 ±0.0004 


-47.61 


1100 




10 


1.1723 ± 0.0005 


-49.05 


1300 




-10 


1.1803 ± 0.0005 


-50.24 


1800 




-30 


1.1899 ±0.0004 


-51.49 


2500 



s 



260 



100 



0.1 



-200 



35 1.1060 ±0.0004 

10 1.1145 ±0.0006 

-10 1.1201 ±0.0005 

-30 1.1224 ±0.0005 

35 1.0514 ± 0.0004 

10 1.0570 ± 0.0005 

-10 1.0555 ± 0.0003 

-30 1.0513 ±0.0006 

35 1.0132 ± 0.0003 

10 1.0158 ±0.0004 

-10 1.0132 ±0.0003 

-30 1.0046 ± 0.0006 

-145 1.022 ±0.009 

35 0.9064 ± 0.0003 

10 0.9212 ±0.0005 

-10 0.9245 ± 0.0003 

-30 0.9283 ± 0.0004 



-47.59 

-49.09 

-50.33 

-51.62 

-47.51 

-49.11 

-50.48 

-51.85 

-47.42 

-49.12 

-50.56 

-52.15 

-56.1 

-46.73 

-48.83 

-50.59 

-52.12 



550 
750 
1500 
2100 
600 
800 
1500 
2100 
600 
800 
1500 
3000 
4000 
600 
800 
1500 
3000 



9 



FIGURES 
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FIG. 1. Structure of liquid water comparing results from experiments [||| and the present simu- 
lations. Each set of curves can be identified as follows (reading from top to bottom at the location 
of the arrow): P = 600 MPa, 465 MPa, 260 MPa, 100 MPa, 0.1 MPa, and -200 MPa. a, The pair 
correlation function h(r) for two of five temperatures studied. Pressures are the same for the ex- 
periments and simulations, with the exception that no experiment was possible at P = —200 MPa. 
Note the pronounced increase in the 3.3 A peak (arrow) when pressure is increased. To facilitate 
comparison with experiments, the simulation temperature is reported relative to the Tmd of the 
SPC/E potential at atmospheric pressure. Similarly, the experimental data are reported relative to 
the Tmd of D2O at atmospheric pressure. Isotopic substitution is reflected by the weighting factors 
used to calculate h(r), given by equation ([[]). Substitution of deuterium for hydrogen is expected 
to have little effect on the individual RDFs. To compare with experimental measurements, we 
use weighting factors for D2O, given by w\ = 0.092, W2 = 0.422, and W3 = 0.486. b, The pair 
correlation function goo( r ) f° r three of five temperatures studied. Note the pronounced increase in 
the 3.3 A peak (arrow) when pressure is increased. Simulations of the ST2 and MCY potential do 
not display a peak at 3.3 A in goo( r ), characteristic of the interpenetrating tetrahedral structure 
expected under high pressure |5],[§. c, The molecular structure factor Sm{q), calculated from the 
Fourier transform of h(r) (Fig. |T|a). Looking at the first peak in Sm{q), the curves are identified 
as described above. Note the shift in the first peak when pressure is increased. 
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FIG. 2. The value qo of the first peak of the structure factor from Fig. |l]c. At low T, both the 
simulated (filled symbols) and experimental (open symbols) jl value of go tend toward the values 
for the two amorphous forms of water, HDA and LDA. 
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FIG. 3. a, Comparison of the structure of the supercooled liquid at AT = —30 °C and 
P = 600 MPa with the experimentally measured structure of HDA solid water from ref. [27]. 
b, Comparison of the supercooled liquid structure at AT = —30 °C and P = —200 MPa with LDA 
solid water pffl. c, Comparison of the structure of the glassy simulation at AT = —145 °C and 



P = 600 MPa with the experimentally measured structure of HDA solid water from ref. [27 \. d, 
Comparison of the glassy simulation at AT = —145 °C and atmospheric pressure with LDA solid 



water [27]. While negative pressure is necessary to observe LDA- like structure in the simulations for 
AT = —30 °C, atmospheric pressure is sufficient to observe LDA-like structure at AT = —145 °C. 
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